# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <errno.h>
# include <string.h>
# include <math.h>
// # include <ctype.h>

extern double kvisc(double T);
extern double dens(double T);
extern double thm_cond(double T);

extern short print_result(unsigned int t_cnt, const struct grid *g)
  {
   int i, j, k, nr, nth, nz;
   char resfile[MAX_FNAME_LEN], err_str[SYSERR_STR_LEN];
   FILE *ptr;
   double ***u, ***v, ***w, *r, *th, *z, _u, _v, _w, v_total, vort_th, div;
   double T_wall, stress_r, stress_th;
#  ifdef SOLVE_ENERGY
   double h, A_total, T_wall_av, h_av, q_av, Nu_av;
#  endif

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'print_result()'...\n");
#  endif

   u = g->u; v = g->v; w = g->w; r = g->r; th = g->th; z = g->z; nr = g->nr; nth = g->nth; nz = g->nz;

/* ============================= Printing result at internal points ========================================== */
   sprintf(resfile, "%02d.%s.%09d.[%E].dat", RUN_CNT, RES_FILE, t_cnt, g->t);

   if((ptr = fopen(resfile, "w")) == (FILE *)0)
     {
      strerror_r(errno, err_str, SYSERR_STR_LEN);
      fprintf(g->tptr, " {'%s' - result file writing error: %s}; ", resfile, err_str); fflush(g->tptr);
      fprintf(g->lptr, "ERROR:: In 'print_result()': Unable to open the result-file '%s' for writing: %s\n",\
                  resfile, err_str);
      return -1;
     }

// Printing header for Tecplot.
   fprintf(ptr, "TITLE = \"Re0 = %E, t = %E Sec", g->Re0, g->t);
#  ifdef SOLVE_ENERGY
   if(g->t < t_HEATER_ON) fprintf(ptr, ", (Heater OFF. It will be ON after %E seconds.)", t_HEATER_ON-g->t);
   else fprintf(ptr, ", (Heater ON since %E seconds)", g->t-t_HEATER_ON);
#  endif

   fprintf(ptr, "\"\nVARIABLES = \"r [mm]\", \"z [mm]\", \"th [degree]\", \"u [mm/s]\", \"v [mm/s]\", \"w [mm/s]\", \"Total Velocity [mm/s]\", \"Vorticity_th [rad/sec]\", \"Divergence [kg/m3.sec]\", \"P [Pa]\"");
#  ifdef SOLVE_ENERGY
   fprintf(ptr, ", \"Temperature[ oC]\"");
#  endif

   fprintf(ptr, "\nZONE I=%d, J=%d, K=%d, DATAPACKING=POINT\n", nr, nz, nth);


   for(j = NthG; j < NthG + g->nth; j++) for(k = NzG; k < NzG + g->nz; k++) for(i = NrG; i < NrG + g->nr; i++)
       {
     /* Evaluating divergence of the velocity field. */
        div = g->density0*DIV_ijk;

     /* Evaluating total-velocity. */
        _u = 0.5*(u[i][j][k]+u[i-1][j][k]);
        _v = 0.5*(v[i][j][k]+v[i][j-1][k]);
        _w = 0.5*(w[i][j][k]+w[i][j][k-1]);
        v_total = pow(_u*_u + _v*_v + _w*_w, 0.5);


     /* Evaluating OMEGA_th, vorticity in theta-direction ('v = 0' is assummed). */
        vort_th = 0.5*( ((u[i][j][k]+u[i-1][j][k])*(z[k+2]-z[k+1])+(u[i][j][k+1]-u[i-1][j][k+1])*(z[k+1]-z[k]))/(z[k+2]-z[k])\
                - ((u[i][j][k-1]+u[i-1][j][k-1])*(z[k+1]-z[k])+(u[i][j][k]+u[i-1][j][k])*(z[k]-z[k-1]))/(z[k+1]-z[k-1]) ) / (z[k+1]-z[k])\
                - 0.5*( ((w[i][j][k]+w[i][j][k-1])*(r[i+1]-r[i+1])+(w[i+1][j][k]+w[i+1][j][k-1])*(r[i+1]-r[i]))/(r[i+2]-r[i])\
                - ((w[i-1][j][k]+w[i-1][j][k-1])*(r[i+1]-r[i])+(w[i][j][k]+w[i][j][k-1])*(r[i]-r[i-1]))/(r[i+1]-r[i-1]) ) / (r[i+1]-r[i]);

        fprintf(ptr, "%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E", 0.5*(r[i]+r[i+1])*1000, 0.5*(z[k]+z[k+1])*1000,\
             0.5*(th[j]+th[j+1])*90/acos(0), _u*1000, _v*1000, _w*1000, v_total*1000, vort_th, div, g->p[i][j][k]);

#       ifdef SOLVE_ENERGY
        fprintf(ptr, "\t%E", g->T[i][j][k]-ABS_ZERO);
#       endif

        fprintf(ptr, "\n");
       }
   fclose(ptr);
/* ----------------------------------------------------------------------------------------------------------------------------- */


   fprintf(g->tptr, "('%s' - written); ", resfile); fflush(g->tptr);


/* ======================================= Printing values at the top impingement wall ========================================= */
   sprintf(resfile, "%02d.top.%s.%09d.[%E].dat", RUN_CNT, RES_FILE, t_cnt, g->t);

   if((ptr = fopen(resfile, "w")) == (FILE *)0)
     {
      strerror_r(errno, err_str, SYSERR_STR_LEN);
      fprintf(g->tptr, " {'%s' - result file writing error: %s}; ", resfile, err_str); fflush(g->tptr);
      fprintf(g->lptr, "ERROR:: In 'print_result()': Unable to open the result-file '%s' for writing: %s\n",\
                 resfile, err_str);
      return -1;
     }

// Printing header for Tecplot.
   fprintf(ptr, "TITLE = \"Re0 = %E, t = %E Sec", g->Re0, g->t);
#  ifdef SOLVE_ENERGY
   if(g->t < t_HEATER_ON) fprintf(ptr, ", (Heater OFF. It will be ON after %E seconds.)", t_HEATER_ON-g->t);
   else fprintf(ptr, ", (Heater ON since %E seconds)", g->t-t_HEATER_ON);
#  endif

   fprintf(ptr, "\"\nVARIABLES = \"r [mm]\", \"theta [degree]\", \"Radial Stress [N/sq.m]\", \"Tangential Stress [N/sq.m]\"");
#  ifdef SOLVE_ENERGY
   fprintf(ptr, ", \"T(top surface) [oC]\", \"Nu\"");
   T_wall_av = q_av = 0; 
#  endif

   fprintf(ptr, "\nZONE I=%d, J=%d, DATAPACKING=POINT\n", nr, nth);

   for(j = NthG; j < NthG + nth; j++) for(i = NrG; i < NrG + nr; i++)
        {
#        ifdef SOLVE_ENERGY
         T_wall = (g->T[i][j][NzG+nz]+g->T[i][j][NzG+nz-1])/2;
#        else
         T_wall = T_INLET;
#        endif

         stress_r = kvisc(T_wall)*dens(T_wall)*(u[i][j][NzG+nz-1]+u[i-1][j][NzG+nz-1])/(z[NzG+nz]+z[NzG+nz-1]);
         stress_th = kvisc(T_wall)*dens(T_wall)*(v[i][j][NzG+nz-1]+v[i][j-1][NzG+nz-1])/(z[NzG+nz]+z[NzG+nz-1]);

         fprintf(ptr, "%E\t%E\t%E\t%E", 0.5*(r[i]+r[i+1])*1000, 0.5*(th[j]+th[j+1])*90/acos(0), stress_r, stress_th);

#        ifdef SOLVE_ENERGY
         h = (g->t < t_HEATER_ON ? 0:(i < NrG+g->nr_heater ? TOP_HEATER_HF/(T_wall - T_INLET):0));

         fprintf(ptr, "\t%E\t%E", T_wall - ABS_ZERO, h*2*r[NrG+g->nr_heater]/thm_cond(T_wall));

         T_wall_av += T_wall*0.5*(r[i+1]*r[i+1]-r[i]*r[i])*(th[j+1]-th[j]);
         q_av += (g->t < t_HEATER_ON ? 0:(i < NrG+g->nr_heater ? TOP_HEATER_HF:0))*0.5*(r[i+1]*r[i+1]-r[i]*r[i])*(th[j+1]-th[j]);
#        endif

         fprintf(ptr, "\n");
        }
   fclose(ptr);
/* ----------------------------------------------------------------------------------------------------------------------------- */

   fprintf(g->tptr, "('%s' - written); ", resfile); fflush(g->tptr);


#  ifdef SOLVE_ENERGY
/* Evaluating average value of "Nu". */
   A_total = 0.5*(r[NrG+nr]*r[NrG+nr]-r[NrG]*r[NrG])*(th[NthG+nth]-th[NthG]);
   T_wall_av /= A_total;
   q_av /= A_total;
   h_av = q_av/(T_wall_av - T_INLET);
   Nu_av = h_av*2*r[NrG+g->nr_heater]/thm_cond(T_wall_av);

   fprintf(g->tptr, "Nu(av) = %E; ", Nu_av); fflush(g->tptr);
#  endif

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'print_result()' ends.\n");
#  endif

   return 1;
  }
